The Figure 2 is the showcase for geomeTriD package to
present multiple genomic signals along with 3D models for a
microcompartment/loop.
library(geomeTriD)
library(geomeTriD.documentation)
library(GenomicRanges)
library(TxDb.Mmusculus.UCSC.mm39.knownGene)
library(org.Mm.eg.db)
library(trackViewer)
library(RColorBrewer)The Region Capture Micro-C data were downloaded from GSE207225. The genomic signals of ChIP-seq were downloaded from GSE178982 and remapped to mm39 genome.
## import Genomic interaction data
geo_acc <- c("DMSO"='GSM6281851',
'IAA'='GSM6281852')
url <- 'https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM6281nnn';
urls <- mapply(geo_acc, names(geo_acc), FUN=function(.f, .cond){
file.path(url, .f, 'suppl',
paste0(.f, '_RCMC_BR1_merged_allCap_', .cond,
'_mm39.merged.50.mcool'),
fsep='/')
})
## import the genomic interactions
range_chr8_wid <- GRanges('chr8:84000000-92000000')
gi <- importGInteractionsFromUrl(urls=urls, resolution=500, range=range_chr8_wid,
format='cool', normalization='balanced')## adding rname 'https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM6281nnn/GSM6281851/suppl/GSM6281851_RCMC_BR1_merged_allCap_DMSO_mm39.merged.50.mcool'## adding rname 'https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM6281nnn/GSM6281852/suppl/GSM6281852_RCMC_BR1_merged_allCap_IAA_mm39.merged.50.mcool'
## create 3d model by mdsPlot function
range_chr8 <- GRanges('chr8:85550000-85800000')
model3d <- lapply(gi, mdsPlot, range = range_chr8, k=3, render = 'granges')## Warning in checkGI(gi, fixedBin = TRUE): There are NA values in the gi score.
## It will be removed.## Some region are missing from the input gi.## initial value 17.000699
## iter 5 value 14.159797
## iter 10 value 13.013113
## iter 15 value 12.575326
## final value 12.445963
## converged## Warning in checkGI(gi, fixedBin = TRUE): There are NA values in the gi score.
## It will be removed.## Some region are missing from the input gi.## initial value 17.114933
## iter 5 value 13.719667
## iter 10 value 12.444197
## iter 15 value 12.102132
## final value 11.986303
## converged
### create gene annotations
#### get all genes
feature.gr <- genes(TxDb.Mmusculus.UCSC.mm39.knownGene)## 2 genes were dropped because they have exons located on both strands of the
## same reference sequence or on more than one reference sequence, so cannot be
## represented by a single genomic range.
## Use 'single.strand.genes.only=FALSE' to get all the genes in a GRangesList
## object, or use suppressMessages() to suppress this message.
#### subset the data by target viewer region
feature.gr <- subsetByOverlaps(feature.gr, range_chr8)
#### assign symbols for each gene
symbols <- mget(feature.gr$gene_id, org.Mm.egSYMBOL, ifnotfound = NA)
feature.gr$label[lengths(symbols) == 1] <-
unlist(symbols[lengths(symbols) == 1])
#### assign colors for each gene
feature.gr$col <- sample(1:7, length(feature.gr), replace = TRUE)
### plot for the detailed region ii to show the open or close loop.
subregion <- GRanges('chr8:85690000-85730000')
model3d.sub <- lapply(model3d, subsetByOverlaps, ranges=subregion)
model3d.sub[['IAA']] <- alignCoor(model3d.sub[['IAA']], model3d.sub[['DMSO']])
### import cohesion ctcf, and yy1 signals, realigned to mm39 for GSE178982
pf <- system.file('extdata', package = 'geomeTriD.documentation')
bws_files <- dir(file.path(pf, 'GSE178982'), '.bw', full.names = TRUE)
(names(bws_files) <- sub('^.*?_(IAA|UT)_(.*?).CPM.*$', '\\2', bws_files))## [1] "CTCF" "RAD21" "SMC1A" "SMC3" "YY1" "CTCF" "RAD21" "SMC1A" "SMC3"
## [10] "YY1"
bw_UT <- bws_files[grepl('_UT_', bws_files)]
bw_IAA <- bws_files[grepl('_IAA_', bws_files)]
signals_UT <- lapply(bw_UT, importScore, format='BigWig', ranges=subregion)
signals_IAA <- lapply(bw_IAA[names(bw_UT)], importScore, format='BigWig', ranges=subregion)
colorSets <- c(CTCF="cyan",YY1="yellow",
RAD21="red", SMC1A="green", SMC3="blue")
for(i in seq_along(signals_UT)){
setTrackStyleParam(signals_UT[[i]], "color", c("gray30", colorSets[names(signals_UT)[i]]))
setTrackStyleParam(signals_IAA[[i]], "color", c("gray30", colorSets[names(signals_UT)[i]]))
}
## set maximal lwd for UT and IAA samples according their max
genomicScoreRanges <- lapply(signals_UT, function(.ele) range(.ele$dat$score))geomeTriD
### create the structure
dmso <- view3dStructure(model3d.sub[['DMSO']], feature.gr = feature.gr,
genomicSigs=signals_UT,
reverseGenomicSigs = FALSE,
genomicScoreRange = genomicScoreRanges,
lwd.maxGenomicSigs = 20,
k = 3, renderer = 'none')## Feature type is missing. Set as gene.
iaa <- view3dStructure(model3d.sub[['IAA']], feature.gr = feature.gr,
genomicSigs=signals_IAA,
reverseGenomicSigs = FALSE,
genomicScoreRange = genomicScoreRanges,
lwd.maxGenomicSigs = 20,
k=3, renderer = 'none')## Feature type is missing. Set as gene.
iaa <- lapply(iaa, function(.ele){
.ele$side <- 'right'
.ele
})
threeJsViewer(dmso, iaa, title = c('DMSO control', 'IAA 3h'))
#widget <- threeJsViewer(dmso, iaa, title = c('DMSO control', 'IAA 3h'), background = 'white')
#tempfile <- 'Fig2.part1.html'#tempfile(fileext = '.html', pattern = 'RCMC_BR1_IAA_vs_DMSO_Klf1_II.3jsViewer.')
#htmlwidgets::saveWidget(widget, file=tempfile)
#utils::browseURL(tempfile)The following plot clearly highlights the top 10 interaction events for DMSO and IAA samples. These interactions shift from the long-distance regions of Hook2 to more localized regions.
## extract backbone coordinates, which will be used as the bone for RCMC data
xyz_dmso <- extractBackbonePositions(dmso)
xyz_iaa <- extractBackbonePositions(iaa)
DMSO <- gi[['DMSO']]
IAA <- gi[['IAA']]
DMSO.subset_II <- subsetByOverlaps(DMSO, subregion)
IAA.subset_II <- subsetByOverlaps(IAA, subregion)
DMSO.subset_II <- DMSO.subset_II[distance(first(DMSO.subset_II), second(DMSO.subset_II))>5000]
IAA.subset_II <- IAA.subset_II[distance(first(IAA.subset_II), second(IAA.subset_II))>5000]
hic_dmso_II <- create3dGenomicSignals(
DMSO.subset_II,
xyz_dmso,
name='dmsoII', # name prefix for the geometry
tag='dmsoII', # name for the layer in the scene
color = c('white', brewer.pal(9, 'YlOrRd')),
topN=10, # only plot the top 10 events ordered by the scores.
lwd.maxGenomicSigs = 3,
alpha=0.5
)
hic_iaa_II <- create3dGenomicSignals(
IAA.subset_II,
xyz_iaa,
name='iaaII', # name prefix for the geometry
tag='iaaII', # name for the layer in the scene
color = c('white', brewer.pal(9, 'YlOrRd')),
topN=10, # only plot the top 10 events ordered by the scores.
lwd.maxGenomicSigs = 3,
alpha=0.5
)
hic_iaa_II <- lapply(hic_iaa_II, function(.ele){
.ele$side <- 'right'
.ele
})
threeJsViewer(dmso, iaa, hic_dmso_II, hic_iaa_II, title = c('DMSO', 'IAA 3h'), background = c('gray10', 'gray20', 'gray20', 'gray10'))
#widget <- threeJsViewer(dmso, iaa, hic_dmso_II, hic_iaa_II, title = c('DMSO', 'IAA 3h'), background = c('gray10', 'gray20', 'gray20', 'gray10'))
#tempfile <- 'Fig2.part2.html'#tempfile(fileext = '.html', pattern = 'RCMC_BR1_IAA_vs_DMSO_Klf1_II.3jsViewer.')
#htmlwidgets::saveWidget(widget, file=tempfile)
#utils::browseURL(tempfile)## R Under development (unstable) (2025-01-08 r87545)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] grid stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] RColorBrewer_1.1-3
## [2] trackViewer_1.43.6
## [3] org.Mm.eg.db_3.20.0
## [4] TxDb.Mmusculus.UCSC.mm39.knownGene_3.20.0
## [5] GenomicFeatures_1.59.1
## [6] AnnotationDbi_1.69.0
## [7] Biobase_2.67.0
## [8] GenomicRanges_1.59.1
## [9] GenomeInfoDb_1.43.2
## [10] IRanges_2.41.2
## [11] S4Vectors_0.45.2
## [12] BiocGenerics_0.53.3
## [13] generics_0.1.3
## [14] geomeTriD.documentation_0.0.1
## [15] geomeTriD_1.1.6
##
## loaded via a namespace (and not attached):
## [1] strawr_0.0.92 rstudioapi_0.17.1
## [3] jsonlite_1.8.9 magrittr_2.0.3
## [5] rmarkdown_2.29 fs_1.6.5
## [7] BiocIO_1.17.1 zlibbioc_1.53.0
## [9] ragg_1.3.3 vctrs_0.6.5
## [11] memoise_2.0.1 Rsamtools_2.23.1
## [13] RCurl_1.98-1.16 base64enc_0.1-3
## [15] htmltools_0.5.8.1 S4Arrays_1.7.1
## [17] progress_1.2.3 plotrix_3.8-4
## [19] curl_6.1.0 Rhdf5lib_1.29.0
## [21] rhdf5_2.51.2 SparseArray_1.7.3
## [23] Formula_1.2-5 sass_0.4.9
## [25] bslib_0.8.0 htmlwidgets_1.6.4
## [27] desc_1.4.3 Gviz_1.51.0
## [29] httr2_1.0.7 cachem_1.1.0
## [31] GenomicAlignments_1.43.0 igraph_2.1.3
## [33] lifecycle_1.0.4 pkgconfig_2.0.3
## [35] Matrix_1.7-1 R6_2.5.1
## [37] fastmap_1.2.0 GenomeInfoDbData_1.2.13
## [39] MatrixGenerics_1.19.1 digest_0.6.37
## [41] colorspace_2.1-1 textshaping_0.4.1
## [43] Hmisc_5.2-2 RSQLite_2.3.9
## [45] filelock_1.0.3 httr_1.4.7
## [47] abind_1.4-8 compiler_4.5.0
## [49] withr_3.0.2 bit64_4.6.0-1
## [51] htmlTable_2.4.3 backports_1.5.0
## [53] BiocParallel_1.41.0 DBI_1.2.3
## [55] biomaRt_2.63.0 MASS_7.3-64
## [57] rappdirs_0.3.3 DelayedArray_0.33.3
## [59] rjson_0.2.23 tools_4.5.0
## [61] foreign_0.8-87 nnet_7.3-20
## [63] glue_1.8.0 restfulr_0.0.15
## [65] InteractionSet_1.35.0 rhdf5filters_1.19.0
## [67] checkmate_2.3.2 cluster_2.1.8
## [69] gtable_0.3.6 BSgenome_1.75.0
## [71] ensembldb_2.31.0 data.table_1.16.4
## [73] hms_1.1.3 xml2_1.3.6
## [75] XVector_0.47.2 pillar_1.10.1
## [77] stringr_1.5.1 dplyr_1.1.4
## [79] BiocFileCache_2.15.0 lattice_0.22-6
## [81] deldir_2.0-4 rtracklayer_1.67.0
## [83] bit_4.5.0.1 biovizBase_1.55.0
## [85] tidyselect_1.2.1 Biostrings_2.75.3
## [87] knitr_1.49 gridExtra_2.3
## [89] ProtGenerics_1.39.2 SummarizedExperiment_1.37.0
## [91] xfun_0.50 matrixStats_1.5.0
## [93] stringi_1.8.4 UCSC.utils_1.3.1
## [95] lazyeval_0.2.2 yaml_2.3.10
## [97] evaluate_1.0.3 codetools_0.2-20
## [99] interp_1.1-6 tibble_3.2.1
## [101] cli_3.6.3 rpart_4.1.24
## [103] systemfonts_1.2.0 munsell_0.5.1
## [105] jquerylib_0.1.4 dichromat_2.0-0.1
## [107] Rcpp_1.0.14 grImport_0.9-7
## [109] dbplyr_2.5.0 png_0.1-8
## [111] XML_3.99-0.18 parallel_4.5.0
## [113] pkgdown_2.1.1 rgl_1.3.16
## [115] ggplot2_3.5.1 blob_1.2.4
## [117] prettyunits_1.2.0 jpeg_0.1-10
## [119] latticeExtra_0.6-30 AnnotationFilter_1.31.0
## [121] bitops_1.0-9 txdbmaker_1.3.1
## [123] VariantAnnotation_1.53.1 scales_1.3.0
## [125] purrr_1.0.2 crayon_1.5.3
## [127] rlang_1.1.5 KEGGREST_1.47.0